library(mvtnorm)

VC <- diag(3)

M <- 1000
n.vec <- c(250, 500, 1000)
alpha.vec <- rbind(c(.1,.1,.05),
                   c(1, 1, .5),
                   c(1.5,1.5,.75))
beta.vec <- rbind(c(1,1,2,0,0,5),
                  c(1,1,2,2,2,1))
                   
                   

set.seed(1234)


for (n in n.vec){
  for (alpha.ind in 1:3){
    for (beta.ind in 1:2){
      
      alpha <- alpha.vec[alpha.ind,]
      beta <- beta.vec[beta.ind,]
      
      for (i in 1:M){
        
        Z <- rmvnorm(n, mean=rep(0, 3), sigma=VC)
        
        eta <- Z[,1]*alpha[1] + Z[,2]*alpha[2] +
          Z[,1]*Z[,2]*alpha[3]
        X <- rbinom(n, 1, pnorm(eta))
        Y <- Z[,2]*beta[1] + Z[,3]*beta[2] + Z[,2]*X*beta[3] +
          (Z[,2]^2)*X*beta[4] + (Z[,3]^2)*X*beta[5] + X*beta[6] + rnorm(n,s=1)
        
        mydata <- data.frame(Z, X, Y)
        colnames(mydata) <- c(paste("z", 1:3, sep=""), "x", "y")
        
        filename <- paste("../data/MCdata.n", n, ".alpha", alpha.ind,
                          ".beta", beta.ind,
                          ".rep", i, ".Rdata", sep="")
        save(mydata, file=filename)
        
      } 
    }
  }
}


## ACE is always 5.00
